Metastability in networks of nonlinear stochastic integrate-and-fire neurons

Neurons in the brain continuously process the barrage of sensory inputs they receive from the environment. A wide array of experimental work has shown that the collective activity of neural populations encodes and processes this constant bombardment of information. How these collective patterns of activity depend on single neuron properties is often unclear. Single-neuron recordings have shown that individual neural responses to inputs are nonlinear, which prevents a straightforward extrapolation from single neuron features to emergent collective states. In this work, we use a field theoretic formulation of a stochastic leaky integrate-and-fire model to study the impact of nonlinear intensity functions on macroscopic network activity. We show that the interplay between nonlinear spike emission and membrane potential resets can i) give rise to metastable transitions between active firing rate states, and ii) can enhance or suppress mean firing rates and membrane potentials in opposite directions.


INTRODUCTION
Populations of neurons can showcase a wide array of complex collective activity patterns that underlie sensory and information processing in the brain.Understanding how these macroscopic patterns of activity emerge from the properties of individual neurons has been a central question in neuroscience.Recently, metastable activitysharp stochastic-or input-driven changes in the firing patterns of the network-has come under investigation as a potential means by which populations could flexibly encode sensory information [1,2].These states are frequently observed as sub-populations of neurons switching between up (high activity) and down (low activity) states [3][4][5][6][7], or transitions between multiple clusters of neurons [8][9][10].Models of up-down transitions commonly describe the activity in these states as bistable attractors in dynamical models of network activity [11][12][13][14][15][16].Similarly, models with multiple clusters can produce multistable attractors that describe metastability between multiple states [17][18][19][20][21].
Classically, theories for describing emergent macroscopic population activity (e.g., the Wilson-Cowan [22,23] or Amari-Grossberg equations [24][25][26]) are commonly understood as coarse grained models of the underlying populations of neurons.These equations are derived under the implicit assumption of a separation of timescales between the synaptic and intrinsic dynamics of neurons [27].Additionally, these models describe the mean population activity, ignoring the role that stochastic fluctuations play in shaping the emergent activity of a population of neurons.Some models introduce stochasticity into firing rate models by introducing an external fluctuating current input, with variance sometimes self-consistently matched to the firing rates of the network to mimic a recurrent Poissonian spiking input.However, this phenomenological addition of noise assumes that correlations in spiking activity are negligible, precluding strong coordination within the network.
To properly understand the influence of recurrent spiking activity, we follow recent modeling efforts that employ stochastic field theory representations of the spiking network dynamics without inserting phenomenological noise to mimic the effects of spiking variability [28][29][30][31][32].In particular, we use a stochastic spiking network that incorporates a hard reset of the neurons' membrane potentials after they emit a spike [33].This hard reset aids in stabilizing network activity, allowing the stochastic neurons to fire guaranteed spikes if their input is large enough, yet remain stable because their membrane potentials are always reset after a spike.
The stability of a network and the possible patterns of collective activity it admits is also intrinsically tied to nonlinearities in neural activity, even at the single neuron level.For example, neurons are often characterized by an intensity (or transfer) function that characterizes their instantaneous rate of firing given their current membrane potential [28,[34][35][36][37][38][39][40].This intensity function cannot be measured directly; instead, the firing rate as a function of mean membrane potential has been estimated by fitting single neuron data.The measured nonlinearities are assumed to be a noisy version of the underlying intensity function, giving insight into the appropriate families of intensity functions to use in theoretical and computational models.For example, in the primary visual cortex the nonlinearities are fit well by threshold-power law functions of the form ⌊V − θ⌋ α + , where ⌊x⌋ + = x if x > 0 and 0 otherwise, and α is an exponent [37][38][39][40].Observed exponents α typically fall in the range of 2 ∼ 5, larger than the rectified linear units often used to model neural firing in mathematically tractable models of stochastic spiking neurons [33,41,42].
The rectifying nature of the threshold power-law stands in contrast with the exponential nonlinearity, an alternative modeling choice that has been used exten-sively in fitting point process generalized linear models (GLMs) and generalized integrate-and-fire models to spiking activity data [29,31,[43][44][45][46][47][48][49][50].The use of exponential nonlinearities in fitting GLMs is primarily for technical convenience, but it motivates us to also investigate whether the collective activity generated by such networks approximates the activity of threshold-power law networks.As shown in Fig. 1, we find that both power law and exponential nonlinearities can fit intensity functions well across cell types and layers of visual cortex in mouse data from the Allen Institute for Brain Science [47], but we will show that the two classes of nonlinearities exhibit several key differences in the collective activity they can generate.
To show this, we use a stochastic field theory formalization [33,[51][52][53][54][55][56] of a biologically realistic model of stochastic leaky and fire (sLIF) neurons with nonlinear intensity functions.We study its dynamics to compare the impact of threshold power-law and exponential intensity functions on population activity.We first calculate the mean-field phase diagram of the stochastic spiking network as a function of its synaptic connectivity strength and external input.We show that superlinear threshold-power law intensity functions, coupled with a hard reset of the membrane potential after spiking, enable metastability between i) a quiescent and active state for subthreshold currents and ii) two active states for super-threshold current.In the latter case the networks can stochastically switch between these two active states.Moreover, we show that the quiescent steady state is a result of the rectifying nature of the intensity.Networks with exponential intensity functions exhibit a monostable regime of continuously varying firing rate or a metastable regime between low and high firing rates.
To capture the impact of fluctuations on population activity, we go beyond mean-field theory to compute fluctuation corrections.The fluctuation corrections to the mean-field theory come from the nonlinear firing intensity and the spike reset.While the spike reset always suppresses activity, the nonlinear firing intensity promotes activity when the nonlinearity is superlinear and suppresses activity when it is sublinear [28].We show that this interplay between the nonlinear intensity and spike reset fluctuations shapes the overall dynamics of the network.
The paper is organized as follows: in Section II A we introduce the model and the corresponding mean-field equations using which we can obtain the phase diagram for the network.In Section II B, II C we obtain the meanfield phase diagrams for the two nonlinear intensity functions and find that they can have at most two stable steady-state solutions.Finally, we analyze the impact of fluctuations on the mean firing rate and membrane potential of the network.This helps us dissociate the impact of the two sources of fluctuation corrections: the spike reset and nonlinear firing intensity.

A. Network model and mean-field equations
We model a network of N stochastic leaky-integrate and fire neurons (sLIF).The membrane dynamics are governed by the stochastic differential equation: where V i (t) is the membrane potential of neuron i at time t, J ij is the connection strength from neuron j to neuron i (equal to J/(pN ) with probability p, the connection probability, and 0 otherwise), and E i is the net external current that each neuron receives.We nondimensionalize the model by setting τ = 1.The spike trains ṅi (t) are conditionally Poisson, with instantaneous firing rate ϕ(V i (t)); i.e., the instantaneous firing rate is a function of the neuron's membrane potential.After the emission of a spike, the last term (− ṅi (t + )V i (t)) resets the membrane voltage to 0; see Fig. 1a).We aim to characterize the steady-states for rectified power-law intensity functions, ϕ(V ) ∝ ⌊V − θ⌋ α + with exponent α > 0 and exponential intensity functions ϕ(V ) ∝ e V −θ , where ⌊x⌋ + = x if x > 0 and 0 otherwise.Here, θ is the activation threshold in power-law networks and a "soft" threshold in exponential networks.In both cases it is not the threshold at which spikes are guaranteed to occur, but above which the probability of spiking increases sharply.
To motivate the use of these particular nonlinear intensity functions, we fit the average membrane potential vs firing rate to individual cells from the openly available electrophysiology data from the Allen Institute for Brain Science [47] (Appendix A).The dynamics of individual cells differs from the model dynamics where the details of the spiking and hyper-polarization are replaced by the hard-reset (Fig. 1b).The fitted exponents fall in the range of 2 ∼ 4 as has been observed previously, with the mean exponent of 3.2 (Fig. 1c-e).
The stochastic dynamics in (1) gives rise to a probability distribution for the membrane potential and firing rate for the network.The joint moment generating functional (MGF) for membrane potentials and spike trains can be expressed as a path integral in the Martin-Siggia-Rose-De Dominicis-Janssen (MSRDJ) formalism, described by the action S[ Ṽ , V, ñ, ṅ]: where x • y = i dt xi (t)y i (t), Ṽ , ñ are called the response variables, which can be interpreted as noise variables driving the membrane potential V and firing rate ṅ, and { j, j, h, h} are "source" fields.Moments of the joint probability distribution can be obtained by taking derivatives of the MGF with respect to the source fields.Generating functionals are a powerful tool for studying stochastic dynamics [54][55][56][57].In neuroscience, these functional methods have previously been applied to firing rate networks, e.g., [42,[58][59][60] and coarse-grained networks of neural activity, e.g., [61][62][63][64][65]. Recent work has applied such methods to spiking networks with "soft resets" that reduce the membrane potential by fixed amounts after a spike [28][29][30].In this work we follow [33] in applying this formulation to networks with hard resets.
For random synaptic connections J ij we can average the MGF over realizations of the synaptic connections to derive an effective dynamical description of the population dynamics.In the weak coupling regime in which the weights J ij scale as 1/N , one can ignore the higher cumulants of the connectivity matrix and describe the collective dynamics of the network using only the mean connection strength J.If we assume that there are M clusters of neurons that are homogeneous within a cluster, the mean-field treatment reduces the N -dimensional dynamics to an effective M -dimensional description of the population statistics within each cluster corresponding to the action where the Greek indices label the different clusters and ⟨ ṅν ⟩ is the population averaged activity that needs to be determined self-consistently from the mean-field equations of motion.This decouples the neurons in the network, such that the average behavior of an individual neuron in a cluster depends only on the mean synaptic field M ν=1 J µν ⟨ ṅν ⟩.This description is exact in the N → ∞ limit [66].The MGF in (3) can then be used to obtain the deterministic mean-field approximation for the population membrane potential and firing rate: Unlike standard dynamic mean-field treatments of firing rate networks, Eq. 4 is not exact.The stochasticity of spiking remains in Eq. 3 and influences the dynamics of the population-averaged membrane potentials V µ (t) through the spike reset ṅµ (t + )V µ (t).Nevertheless, deterministic mean-field approximations often paint a good qualitative picture of the dynamics of the system, and we will first investigate the different phases of network activity within this approximation.Then, in Section II D, we account for the influence of spiking fluctuations on network activity by expanding around steady-state solutions of the mean-field equations to calculate Gaussian (one-loop) corrections.Alternatively, we can take advantage of the hard reset to estimate the exact firing rate of the networks using renewal theory.For self-averaging networks, the density in (3) factorizes over clusters.Each neuron (µ ∈ M ) spikes independently given a self-consistent mean-field input.When the network reaches a steady-state, the selfconsistent input becomes independent of time, in which case the spiking dynamics can be treated as a renewal process due to the hard reset of the membrane potential after spike emission.After each spike, the membrane potential asymptotically evolves to the net input from the reset potential (=0 here): where is the net input that the neuron receives with ⟨ ṅν ⟩ to be determined selfconsistently.For an arbitrary nonlinearity ϕ(V (t)) the inter-spike interval density p µ (s) is the product of the instantaneous firing probability at time s and the survival probability until time s: The mean inter-spike interval for each cluster is ⟨s µ ⟩ = ∞ 0 ds s p µ (s).The self-consistent firing rate is the inverse of the mean inter-spike interval: ⟨ ṅµ ⟩ = 1/⟨s µ ⟩, and the roots of this system of equations allow us to estimate the steady-state firing rates for the network.

B. Homogeneous Networks
We're interested in examining the steady states of the mean-field equations (4).For simplicity, we start with a homogeneous excitatory network of N coupled neurons (Fig. 2a).In this case, the mean-field theory is onedimensional and a phase diagram can be obtained by investigating the fixed points of the mean-field dynamics.The phase diagram will be a function of the external input (E) and the mean coupling strength (J).

Threshold Power-Law Intensity
We first investigate population activity of threshold power-law networks.The mean-field equation of motion simplifies to where we have fixed the threshold at θ = 1 without loss of generality.The intensity function is zero below threshold so V = E is the only possible sub-threshold solution, which does not exist if E > 1. Above threshold, one can show that for any value of the exponent α, the above equation can have a maximum of 3 steady-state solutions, giving rise to bistability in the mean-field dynamics (see Appendix B).Bistability in the mean-field dynamics corresponds to metastability in the stochastic spiking network dynamics.It was previously shown in Ref [33], that with a threshold linear intensity function (α = 1), the homogeneous network can either be monostable or metastable (Fig. 2b).It is not possible to explicitly calculate the full phase diagram solution for arbitrary α, but the cases α = 2, 3 are analytically tractable and we give the full solutions here.For α > 3 we can still characterize several general properties of the phase diagram.
We begin with the case α = 2.The mean-field equation of motion ( 4) for the threshold-quadratic intensity function is a cubic equation (assuming V is above threshold): While one can solve for the roots of this equation, we are more interested in the phase diagram of the network in the E−J plane, which motivates us to find the boundaries between cases for which we have multiple real solutions above threshold.We can identify the phase boundaries by analyzing the discriminant, which determines the boundary between a single real solution and three real solutions of Eq. ( 8).The discriminant is a quadratic in E whose roots are functions of J, yielding the mean-field boundaries between bistable and monostable phases; see Appendix B for details.Analogous to the threshold-linear intensity function, the mean-field dynamics are bistable for sufficiently large values of J.However, in threshold-linear networks the bistability is only between a quiescent state and an active state that exists only for subthreshold input E. The superlinear networks retain this regime, but can also exhibit bistability for superthreshold input, and in this region, the bistable state has two active steady-state solutions (Fig. 2c).This allows for stochastic switching between two corresponding metastable active states in the full stochastic network (Fig. 2f).
Turning to α = 3, the mean-field equation of motion (4) for the threshold-cubic intensity function is a quartic polynomial when V is above threshold: where the coefficients {b, c, d, f As in the previous case, we are interested in identifying phase boundaries, which we can again identify by investigating the discriminant of the quartic equation, this time a cubic in E; we give the details in Appendix B. We again find that only two stable fixed points exist above threshold, meaning the bistable state can have two active steady-state solutions for superthreshold input (Fig. 2d).Furthermore, we observe that increasing the exponent α from 2 to 3 expands the region in which both the bistable states are active (Fig. 2c vs d).
For the homogeneous network, the membrane dynamics can be treated as a renewal process where after each spike the membrane potential evolves as: The mean firing rate can then be estimated by numerically finding the roots of (see Appendix C for details): where and 2 F 1 (a, b, c, x) is the hypergeometric function.(The subscript p refers to the power-law intensity function.)Fig. 2g shows the firing rate predictions of the meanfield theory, the one-loop correction (discussed in Section II D), and renewal theory compared to simulations.
While analytic results are difficult to obtain for α > 3, we can show that there are at most two stable fixed points above threshold, meaning the mean-field dynamics can be at most bistable (Appendix B).We can also determine the coordinates of the 'cusp' of the phase diagram for arbitrary α, finding that it traces a non-monotonic path through the J−E plane (Fig. 2h).Our result suggests that for α ≳ 2 an increase in exponent of the intensity function expands the region of bistability where two active metastable steady-states can coexist.
In all cases the bistable region of the (J, E) plane is subdivided in two.Below the threshold E = 1, there exists a quiescent steady state and a high-rate steady state.Above threshold, the quiescent state becomes a low-rate state.Both active states are metastable in the full stochastic network, which can stochastically transition between them (Fig. 2f).If the states are far apart a transition can be induced with the application of an external current (Fig. 2e).

Exponential intensity functions
For an exponential intensity function, ϕ(V ) = e V −θ , we find that the mean-field dynamics can be bistable between two states or monostable, similar to the threshold power law networks.However, the exponential nonlinearity is never zero for any finite input so the network does not admit a quiescent state.Instead, in the monostable regime the firing rate changes continuously as the parameters J or E are tuned.We show in Appendix D that the mean-field dynamics of the homogeneous exponential network are bistable if: where W 0,−1 ≡ W 0,−1 (−e θ+1−J ) is the Lambert W function. Figure 3a shows the mean-field bistability boundaries for various values of θ ∈ [0, 5].We see these states in the stochastic spiking network as well, where the network can either be monostable (Fig. 3b) or metastable (Fig. 3c).
We can also estimate the exact firing rates for the exponential intensity using renewal theory.In this case, the self-consistent firing rates are roots of (see Appendix C for details): where x dt e −t /t is the exponential integral.Figure 3d shows the firing rate predictions for a homogeneous network with an exponential intensity function.

C. Excitatory-Inhibitory Networks
Neuronal networks in the brain are not homogeneous.Realistic networks follow Dale's law: there are both excitatory and inhibitory neurons in these networks.To look at more biologically realistic networks, here we consider an excitatory-inhibitory (EI) network with pulse coupling (Fig. 4a).We introduce an additional parameter (g) that quantifies the relative strength of inhibition such that the mean connectivity matrix has the following structure: If both populations receive the same input, the mean-field theory for this EI network reduces to the one-dimensional model with the replacement J → J(1 − g) in the homogeneous network results.
In Fig. 5 we show the results for threshold power law intensities in the E−g plane (Fig. 4b), along with a raster plot of a simulation to confirm the full stochastic network exhibits metastability (Fig. 4c).We compare the mean-field predictions of the firing rates with the oneloop corrections (discussed in the next section), and the estimates from renewal theory and the simulations (Fig. 4d).
In Fig. 6 we show the corresponding results for E-I networks with exponential nonlinearity, this time in the J−g plane (Fig. 5b), along with demonstration of inputdriven metastability (Fig. 5c) and the comparisons of mean-field, one-loop, renewal theory, and simulation estimates of the firing rates (Fig. 5d).

D. Impact of fluctuations on the mean activity
While the mean-field theory provides a reasonably accurate qualitative description of the network dynamics and phase transitions in the network, it is quantitatively inaccurate (e.g., Figs.2g, 3-5d) because it ignores all fluctuations.These fluctuations impact the mean activity due to i) the nonlinear spike reset and ii) the nonlinear intensity function.In isolation, these have opposite effects: the spike reset suppresses activity [33], while a concave-up intensity function promotes activity [28].We next ask: in the presence of both sources of nonlinearities, do fluctuations enhance firing or suppress it?
To answer this we calculate the so-called perturbative 1-loop corrections to the mean-field predictions [33].In the path integral formalism these corrections can be evaluated by perturbatively expanding the MGF under a Gaussian approximation.This can be achieved by the use of a diagrammatic expansion (Appendix E).These Feynman diagrams are constructed from edges and vertices.
The edges in the Feynman diagrams correspond to the linear responses of configuration variables (V, ṅ) to the response variables ( Ṽ , ñ) -also called propagators.Since we have 2 each of configuration and response variables, they give rise to 4 propagators in the model, given here in the frequency domain: where ϕ (n) (V ) is the nth derivative of the intensity and V and n are the mean-field estimates of the membrane potential and firing rate, respectively (Eq. 4 for a single population).
The vertices in the Feynman diagrams correspond to source and interaction terms in the action (Appendix E).The source terms correspond to sources of fluctuationshere, the stochastic spike emission.The interaction terms are generated by the nonlinearities, i.e., the spike reset and the nonlinear intensity function.The vertices that give rise to 1-loop corrections are: where the first vertex is a source vertex that emits 2 response variables ñ and has an amplitude ϕ (0) ( V ), the second vertex comes from the nonlinear spike reset and has an amplitude −1, and the last vertex comes from the nonlinear intensity function and has an amplitude ϕ (2) ( V )/2.The second and third vertices are interaction vertices.Using the Feynman rules derived by [33], we can evaluate the loop corrections to V and ṅ.Both quantities get corrections from the intensity function as well as the spike reset.The perturbative 1-loop corrections to the mean-field theory are: which can be evaluated using the residue theorem (see [30,56,67] for an introduction to evaluating frequency integrals): As noted earlier, the corrections due to fluctuations in (18), (19) come from two terms.The first correction comes due to the spike reset, and always suppresses both firing rate and membrane potential.The second correction is due to the nonlinear intensity function.It can either suppress or promote the firing rate/membrane potential depending on the concavity of ϕ.Furthermore, the second correction appears with opposite sign in the two equations.This is because the linear response of the rate to spike fluctuations is is positive, while the linear response of the voltage to spike fluctuations is negativeanother manifestation of the spike reset.That term could suppress firing rates while paradoxically increasing the mean membrane potential, or vice versa.This raises the question: under what conditions do fluctuations promote or suppress the mean membrane potential and firing rate of the network?
We find that fluctuations can promote the membrane potential only when the intensity function is sublinear and promote the firing rate only when the intensity function is superlinear.The boundary between enhancement versus suppression occurs when the two fluctuation corrections in ( 18) and ( 19) cancel.With α > 1 there are two mutually exclusive possibilities: the corrections to the mean membrane potential vanish or the corrections to the firing rates vanish.These occur when: These curves define the surfaces where the one-loop fluctuation corrections vanish (Fig. 6a).These can be viewed in the E−J plane by substituting the above solutions for V in the mean-field equations of motion ( 4): ) For the monostable regimes, these are lines in the E−J plane that separate the regions where the membrane potential/firing rate are promoted/suppressed.For the sublinear intensity functions, increasing α increases the intercept of the lines (Fig. 6b), while for superlinear intensity functions, increasing α increases the magnitude of the slope as well as intercept of the lines (Fig. 6c).The equations suggest that for a stronger coupling and a larger external input, fluctuations always suppress the firing rate and the membrane potential.This is a result of the nonlinear spike reset that prevents the dynamics from exploding despite the unbounded nonlinearity.Fig. 6d shows the difference between the firing rate of the stochastic spiking network and the mean-field prediction of the firing rate.
For the EI network, the same boundaries can be estimated with the replacement J → J(1 − g) in (21).The boundaries for superlinear intensity functions and the comparisons with simulation are shown in Figure 6e,f.The equations suggest that-similar to the homogeneous case-fluctuations always suppress the firing rates and membrane potentials for the parameter regimes where the network activity becomes very large.

Fluctuation corrections for exponential intensities
The exponential intensity functions receive the same fluctuations corrections as in (18) and (19).For an exponential, however, ϕ (2) (V ) > 0 for all V > −∞.Hence, there are no parameter regimes where the membrane potential can be promoted by the fluctuations, since both correction terms to V come with a negative sign.As before, we assume the intensity function has the form ϕ(V ) = e V −θ where θ is the soft threshold of the intensity function.The nullcline where the fluctuations corrections vanish for an exponential intensity function is: which again is a straight line in the E−J plane (Fig. 7a).
For the EI network, the same boundaries can be estimated with the replacement of J → J(1−g) in ( 22) (Fig. 7b).Fig. 7c, d show the difference between the firing rates of the stochastic spiking network and the meanfield prediction for the homogeneous population (E−J plane) and the EI populations (J−g plane) for fixed E, respectively.

III. DISCUSSION
We analyzed how single neuron intensity functions impact the overall network activity in a network of stochastic leaky-integrate and fire neurons.Here, we specifically focused on the threshold-power law intensity functions, which have been fit well to experimentally recorded activity in the visual cortex [37][38][39][40], and the exponential intensity functions that have been extensively used in fitting point process GLMs to neuronal recordings [29,31,[43][44][45][46][47][48][49][50].We demonstrated that networks with superlinear threshold-power law intensity functions can have a maximum of 2 stable steady states in the weak coupling regime (J ij ∼ 1/N ).Similar to threshold linear intensity functions [33], the homogeneous and EI networks can be monostable (active or quiescent), or metastable between a high firing rate and quiescent state.However, for superlinear intensities a new metastable region that admits active high-and low-firing rates emerges [11,68], reminiscent of the up-down transitions seen in numerous cortical areas [3][4][5][6][7].We find that increasing the exponent of the intensity function continuously expands the region where these active states co-exist.(Fig. 2).This high-low metastability also appears in networks with exponential intensity functions (Fig. 3), owing to the fact that the firing intensity is never zero-the "subthreshold" firing probability is small but still finite.Depending on the network parameters, the transitions between the two active states can be driven stochastically, or with the application of an external current.
While the mean-field theory provides a decent qualitative description of the network dynamics, it neglects all fluctuations.Using a diagrammatic expansion, one can add the effects of nonlinear fluctuations on the network activity.In the sLIF model, the corrections to the mean-field come from the nonlinear intensity function and the hard reset of the membrane potential after a spike [33].The effect of the nonlinear intensity depends on its concavity.For the threshold power law intensity, when the intensity function is sublinear (α < 1) fluctuations always suppress the firing rate and promote the membrane potential.When the intensity function is superlinear (α > 1)-thought to be the more relevant case experimentally-fluctuations always suppress the firing rate and promote the membrane potential.The two corrections to both variables (firing rate and membrane potential) can vanish simultaneously only when the intensity is threshold linear α = 1.For the exponential intensity, the concavity is always positive, so that the nonlinear intensity always suppresses the membrane potential and always promotes the firing rate.Stochastic transitions in multistable firing states have been ubiquitously observed in experimental recordings, and are often associated with behavioral readouts in the tasks [3-5, 9, 10, 12, 16, 18-20].Understanding how networks parameters can shape such stochastic transitions in multistable networks of spiking neurons has been a long standing question in neuroscience [8,17].The field theo-retic formulation provides a powerful method to investigate these transitions.Previous work has estimated transition probabilities in networks of binary units (loosely interpreted as active or inactive neurons) [21,62,69].Adaptations of this formalism can be potentially be applied to these networks of spiking neurons to delineate the underlying mechanisms that drive these stochastic transitions and elucidate how network parameters shape the transition rates between different states.We leave this as a direction for future work.
Another avenue for extension of our work is the investigation of networks with "strong" synaptic weights that scale as 1/ √ N rather than 1/N .In strongly coupled enetworks the variance of the inputs to a neuron is O(1), and typically balanced in such a way that the mean input cancels out.Firing rate networks with strong synaptic connections have been studied extensively [42,58,60,70,71], but extending these methods to networks of spiking neurons is non-trivial due to the feedback between the spiking and membrane potential evolution.
Extending the analysis to a finite-network size and strong synaptic connections would yield valuable insights about the mechanisms that drive stochastic transitions in networks of spiking neurons.where γ is the gain of the intensity.The details for the fit are as follows: 1.For each level of stimulation, we extracted the mean membrane potential of the cell, excluding the membrane potential values between the time of reaching the threshold (extracted from the spike feature data) and the time of reaching the trough (trough time).
2. We estimated the firing rate by counting the number of spikes during the stimulation divided by the total stimulation time.
3. For the power law intensity, we fixed the threshold at the maximum value of the mean membrane potential where the firing rate was zero.We then fit the gain (γ) and exponent (α).
4. For the exponential intensity, we treated both gain (γ) and threshold (θ) as free parameters in our fit.
Appendix B: Steady-states for Threshold Power Law Intensity Function Let the right-hand-side of the mean-field equation of motion, assuming a superthreshold mean membrane potential, be where we have set the threshold to θ = 1.The values of V where f ( V ) vanishes correspond to the steady states of the mean-field dynamics.The solutions cannot generally be obtained in closed form, but we can show that at most two solutions exist that satisfy V ≥ 1.We do so by investigating the points of inflection of f ( V ), looking at the zeros of the second derivative: ) From this we see that there are only two points of inflection: Vinf,1 = 1 or Vinf,2 = (α−1)J+2 α+1 .Thus, there is only one point of inflection above threshold for a general α.We also see that lim V →∞ f (V ) → −∞.If E > 1, the value of the function at Vinf,1 , f ( Vinf,1 ) is positive, which implies the function can have a maximum of 3 roots above threshold.If E < 1 then f ( Vinf,1 ) < 0, and the function can have a maximum of 2 roots above threshold.In both cases ones of those roots is an unstable fixed point of the mean-field dynamics, and will not be observed in practice.Accordingly, the mean-field steady-state equation can have a maximum of two stable steady states.
We can also use this to determine the location of the critical point where the two roots meet for a general α.The critical point is the point where the first three derivatives of the steady-state equation vanish simultaneously From this we see that lim α→1 [J, E] = [2.0,1.0], and lim α→∞ [J, E] = [2.0,2.0].The critical point continuously increases for α > 1, thus, expanding the region above threshold where the active-active states can coexist.
For specific values of α we can obtain exact solutions for the mean-field phase boundaries.In the case α = 2 the roots of the fixed point equation can be shown to have 1 or 3 real-valued roots.The transition from 1 to 3 can be determined by investigating the discriminant, a quadratic equation in E: where A(J) = [32(J + 2)(J + 1) − 4(J + 2) 3 − 54J] and B(J) = [36J(J + 1)(J + 2) − 4J(J + 2) 3 + 4(J + 2) 2 (J + 1) 2 − 32(J + 1) 3 − 27J 2 .The roots of this quadratic equation give us the mean-field boundaries that separate bistable and monostable phases.For α = 3 the root equation is quartic in V , where the coefficients {b, c, d, f } are b = −(J + 3), c = 3(J + 1), d = −3J, f = (J + E).The discriminant is cubic in E: where . The roots of this cubic equation, along with the condition that V must be above threshold and the following inequalities: define the boundaries separating the bistable and the monostable phases.
Appendix C: Exact firing rates from renewal theory After averaging over the synaptic weights J ij and neglecting variance in the weights, which is negligible when N is large, we obtain a population-averaged stochastic differential equation that follows from the action Eq. ( 3): where ⟨ ṅµ (t)⟩ is the population-averaged input to each cluster.In the large N limit this concentrates to a constant value that must be determined self-consistently.In this limit the clusters are formally decoupled and interact only through the self-consistent input.Due to the hard reset of the membrane potential after a spike, each cluster's spike train can be treated as a renewal process.Between spikes the membrane dynamics of each cluster evolves independently according to: where C µ = E µ + M ν=1 J µν ⟨ ṅν ⟩ is the net input that drives the membrane potential of the cluster µ.The instantaneous firing intensity, from the model definition is simply ϕ(V µ (t)).The inter-spike interval density is then the product of the instantaneous firing intensity and the survival probability independently for each cluster: The mean inter-spike interval can then be evaluated as ⟨s µ ⟩ = ∞ 0 ds sp µ (s), and the mean firing rate is the inverse of the mean inter-spike interval ⟨ ṅµ ⟩ = 1/⟨s µ ⟩.For the nonlinearities we consider here, it is not possible to get a closed form expression for the density.However, we can take a semi-analytic approach to estimate the exact firing rates for the two nonlinearities using renewal theory.In the following we drop the µ subscripts in the intermediate results for simplicity, as different clusters interact only through the coefficients C.

Threshold Power Law Intensity
The inter-spike interval density for the threshold power law intensity function is: While it is not possible to get a closed form expression for this integral, the integral in the exponential of Eqn.(C5) can be evaluated in closed form: To get a closed form expression for this integral, we make a change of variables to y = C(1 − e −t ) − θ, after which the integral can be identified as proportional to a representation of the hypergeometric function: The mean firing rate can then be estimated by numerically finding the roots of the following system of equations, restoring the cluster subscripts and recalling that . (C7)

Exponential Intensity
The inter-spike interval density for the exponential intensity function is: Again, the integral in the exponential of Eqn.(C8) can be evaluated in closed form by making a change of variables y = Ce −t , identifying the result in terms of the "Exponential Integral" special function: ) where Ei(x) = − ∞ −x dt e −t /t is the exponential integral.The mean firing rate can then be estimated by numerically finding the roots of the following system of equations (after restoring the explicit subscipts): Appendix D: Steady-states for exponential intensity functions The steady-state equation for the exponential intensity function is: We can find the extremum of the function by demanding that its derivative vanishes: where, W 0,−1 are the two real branches of the Lambert W function.Thus, the derivative can vanish at either 2 points, or does not vanish at all.Since the argument is the negative of an exponential, the derivative can vanish at exactly 2 points if −1/e < −e θ+1−J < 0; otherwise the function is strictly decreasing.Therefore, the necessary condition for the derivative to vanish at 2 points is: Since the function can have either 2 extrema or no extrema at all, the mean-field equations can have 3 real roots if the values of the function f (V ) at V0,−1 have opposite signs, and will have only 1 real root if the values of the function f (V ) at V0,−1 have the same sign or if f ′ ( V ) never vanishes (J < θ + 2).Substituting the solution (D3) in (D1): f ( V0,−1 ) = (E − J) + 1 − W 0,−1 1 + e J−1−θ+W0,−1 (D5) where we have suppressed the argument of Lambert W function.Thus, there exist 3 steady-states in the network if f ( V0 ) > 0 and f ( V−1 ) < 0, or: The mean-field dynamics are bistable between the two curves above, and monostable otherwise.

Appendix E: Feynman rules for fluctuation corrections
If the fluctuations are weak, we can expand the MGF in (3) in a perturbative series that adds corrections to the mean-field predictions.This can be formally achieved by splitting the action into two parts: i) the free part (S F ) with terms that are up to bi-linear in the fields [ Ṽ , V, ñ, ṅ]; ii) the interaction part (S I ) containing the terms with higher powers of the corresponding fields.Splitting the fields in a background mean and fluctuations (V → V + δV, ṅ → n + δ ṅ), we can find the free and interacting part of the action in the MGF for the homogeneous case in (3): The MGF can now be expanded around the solution to the mean-field theory in the first line of the free action above: The perturbative expansion uses the fact that the case of a Gaussian action (S F ) can be solved exactly.The interacting action (S I ) is then expanded in a functional Taylor series around this Gaussian solution to systematically add the effect of non-Gaussian features of the distribution.An arbitrary moment can be evaluated using: where Z[0] is the MGF evaluated at vanishing source terms, equal to 1 because of normalization of the probability density.At the lowest order in the perturbation, the expectation value of the fluctuations is zero (⟨δV ⟩ = ⟨δ ṅ⟩ = 0), so that (⟨V ⟩ = V , ⟨ ṅ⟩ = n).This is the mean-field solution.The loop expansion gives these fluctuations a non-zero expectation which provides fluctuation corrections to the mean-field solution.The MGF can now be expanded around the free part of the action to calculate expectations beyond the lowest order:
(≥ 1, 0) (0, ≥ q) TABLE II.The vertices in Feynman diagrams correspond to the terms in the interacting part of the action (SI ).Each vertex can have incoming field variables (in-degree) or outgoing response variables (out-degree).
Propagator Edge Factor (t, s) ∆nñ(t, s) δ(t − s) − V ϕ (1) e −(1+n+ V ϕ (1) )(t−s) Θ(t − s) ∆ n Ṽ (t, s) ϕ (1) e −(1+n+ V ϕ (1) )(t−s) Θ(t − s) ∆ V Ṽ (t, s) e −(1+n+ V ϕ (1) )(t−s) Θ(t − s) ∆V ñ(t, s) − V e −(1+n+ V ϕ (1) )(t−s) Θ(t − s) This expansion can be organized diagrammatically using Feynman rules to estimate various moments of a given distribution in terms of the propagators evaluated at the mean-field solution.Using Wick's theorem, any expectation can be broken into a sum over product of propagators ⟨x(t)x(t ′ )⟩ → [∆ nñ , ∆ V ñ, ∆ n Ṽ , ∆ V Ṽ ] (where x = V or ṅ and x = Ṽ or ñ).The Ito condition ⟨x(t)x(t)⟩ = 0 ensures that the propagators are purely causal in the time domain.Thus, the only terms that survive in the above expansion are the terms that have an even number of response [ Ṽ , ñ] and field [V, ṅ] variables, and where each variable in these pairs has a different time argument (comes from different term in the expansion).The diagrams are constructed out of vertices and edges.The vertices (Table II) are the terms in the interacting part of the action (S I ), and the edges (Table III) are constructed out of the propagators evaluated at the solution to the mean-field theory.The diagrammatic rules to evaluate an arbitrary expectation ⟨δ ṅ(t 1 ) . . .δ ṅ(t a )δV (t a+1 ) . . .δV (t a+b )⟩ are as follows: 1. Place a + b external vertices corresponding to the expectation value that needs to be evaluated.
2. Use the vertices and edges in Tables II, III to construct all possible connected graphs such that the external vertices have only one incoming propagator.
3. The analytic expression for the diagram can be obtained by multiplying all the factors of edges and vertices together.In the frequency domain, each edge is assigned its own frequency variable and the expression can be evaluated by integrating over all internal frequencies.In the time domain, each vertex carries its own time argument, and the propagators come with the time arguments of the vertices they connect.The expression can be evaluated by again integrating over all internal time points.
4. The resulting expectation is the sum of all such diagrams that can be constructed in (2).
For the action in (E1), the vertices and the resulting 1loop diagrams are (diagrams generated with [72]): The time arguments corresponding to each vertex are: ( -t 1 ), ( , -t 2 ).The above diagrams in terms of the propagators and edges are: (E7) These expression can be evaluated either in the time domain or can be transformed into the frequency domain and evaluated using the residue theorem.For a detailed introduction to expansions of the moment generating functions see [56,67], path integrals specifically applied to stochastic differential equations see [54], calculations involving loop integrals in the frequency domain see [30], perturbative expansion for this specific model and its connections to the effective action see [33].

FIG. 1 .
FIG. 1. a) Voltage trace showcasing multiple hard-resets after a spike is emitted by a neuron in the model.b) Voltage trace of an example neuron from the cell electrophysiological recording showcasing the full spike generation and post spike hyper-polarization.The model differs from the true dynamics of the cell since the details of the spike generation and hyper-polarization are replaced by the hard-reset.c) Mean membrane potential vs firing rate data for an example cell (empty red circles), with the power law (pink) and exponential (dark blue) fits.The text shows the estimated parameters for both fits.d) The average estimated exponent of the threshold power law fit for different CRE lines.The inset shows the estimated gain for the fit.e) Same as d but for the exponential fit showing the threshold and gain (inset).

FIG. 2 .
FIG. 2. a) Schematic for the homogeneous network.The network has a coupling strength J and receives an external current E. b) Mean-field (MF) phase diagram for the homogeneous network with α = 1 in the E−J plane, separating the bistable (B-H, Q) high firing rate and quiescent region from the monostable high (H) firing rate and quiescent (Q) regions.c) Same as b for α = 2 separating the bistable (B-H, Q) high firing rate and quiescent, bistable (B-H, L) high and low firing rate regions from the monostable high (H) firing rate and quiescent (Q) regions.d) Same as c for α = 3. e) Raster plot for the stochastic spiking network (α = 2) at the parameter marked with cross in panel c (J = 3.2, E = 1.05), illustrating input driven transition between the two active states.f ) Raster plot for the stochastic spiking network (α = 2) at the parameter marked with square in panel c (J = 3.0, E = 1.07), illustrating stochastic transition between the two active states.g) MF (green), 1-loop (cyan) and renewal theory (dark blue) firing rate predictions compared to simulations (brown) for fixed J = 3.0 (α = 2).Inset highlights the low firing rate state in the black square.h) MF phase diagram that depicts that the critical point continuously increases with the exponent of the intensity and reaches a limiting value for arbitrarily large α.Violet: α = 1, Blue: α = 2, Yellow: α = 3, Red: α → ∞.

FIG. 3 .
FIG. 3. a) MF phase diagram for the exponential intensity function separating the bistable (B) regime from the monostable (M) regime for various values of θ ∈ [0, 5] b) Raster plot for the stochastic spiking network (θ = 1) at the parameter marked with square in panel a (J = 4.0, E = −0.75),illustrating the monostable states.c) Raster plot for the stochastic spiking network (θ = 1) at the parameter marked with cross in panel a (J = 4.0, E = −2.0),illustrating input driven transition between the the two active states.d) MF (green), 1-loop (cyan) and renewal theory (dark blue) firing rate predictions compared to simulations (brown) for fixed J = 6.0 (θ = 3).Inset highlights the low firing rate state in the main plot.

FIG. 4 .
FIG. 4. a) Schematic for the excitatory-inhibitory (EI) network (E: Orange, I: Blue).The excitatory cluster has a coupling strength J, the relative strength of inhibition is g, and both populations receive an external current E. b) MF phase diagram for threshold power law intensity in the E−g plane for α = 1 (dashed), α = 2 (dot dashed) and α = 3 (solid).The phase are same as labelled in Fig. 2. c) Raster plot for the stochastic spiking network (α = 2) at J = 5.0, g = 0.4, and E = 1.09 illustrating stochastic transition between the two active states.d) MF (green), 1-loop (cyan) and renewal theory (dark blue) firing rate predictions compared to simulations (brown) for fixed J = 5.0, g = 0.4 (α = 3).Inset highlights the low firing rate state in the main plot.
FIG. 5. a) MF phase diagram for the exponential intensity in the E−g plane for various values of θ ∈ [0, 5].The phase are same as labelled in Fig. 3. b) Same as a but in the J−g plane for θ = 1.c) Raster plot for the stochastic spiking network (θ = 1) at the parameter marked with cross in panel a (J = 8.0, g = 0.4, and E = −2.5),illustrating input driven transition between the the two active states.d) MF (green), 1-loop (cyan) and renewal theory (dark blue) firing rate predictions compared to simulations (brown) for fixed J = 8.0 and g = 0.4 (θ = 1).The inset shows the low firing rate state in the main figure.

FIG. 6 .
FIG. 6. a) Phase diagram in the V −α plane showing the regions where fluctuations promote/suppress the mean firing rate and membrane potential.The curves for α < 1 and α > 1 meet at the same point as α → 1, where the corrections due to the intensity vanish.b Phase diagram in the E−J plane depicting the lines where the fluctuations vanish for sublinear intensity functions (α ∈ [0, 0.9]).Fluctuations promote the membrane potential to the left of each line and suppress the membrane potential to the right of each line.c) Same as b but for superlinear intensity functions (α ∈ [1.5, 5.0]).Fluctuations promote the firing rate to the left of each line and suppress the firing rate to the right of each line.d) The difference between simulation and mean-field predictions of the firing rate for a homogeneous network with a superlinear intensity (α = 2).The solid black line represents the theoretical predictions for the points where fluctuations vanish.e) Same as c but in the J−g plane for an EI network.Here, fluctuations promote the firing rate to the right of each line and suppress the firing rate to the left of each line.f ) Same as d but in the J−g plane for fixed E = 1.5.

FIG. 7 .
FIG. 7. a) Phase diagram in the J−g plane for fixed E = 1.5 depicting the lines where the fluctuations vanish for exponential intensity functions (θ ∈ [0, 5]).Fluctuations promote the firing rate to the left of each line and suppress the firing rate to the right of each line.b) The difference between simulation and mean-field predictions of the firing rate for a homogeneous network with exponential intensity (θ = 1).The solid black line represents the theoretical predictions for the points where fluctuations vanish.c) Same as a but in the J−g plane for for an EI network.Here, fluctuations promote the firing rate to the right of each line and suppress the firing rate to the left of each line.d) Same as b but in the J−g plane for fixed E = −0.5.

TABLE III .
The edges in the Feynman diagrams correspond to the propagators evaluated from the bilinear part of the free action (SF ).Each propagator measures the linear response a field variable to perturbations in the other variable, and are causal due to the Ito condition.Here Θ(t − s) is the Heaviside step function.